Catch bond models may explain how force amplifies TCR signaling and antigen discrimination

The TCR integrates forces in its triggering process upon interaction with pMHC. Force elicits TCR catch-slip bonds with strong pMHCs but slip-only bonds with weak pMHCs. We develop two models and apply them to analyze 55 datasets, demonstrating the models’ ability to quantitatively integrate and classify a broad range of bond behaviors and biological activities. Comparing to a generic two-state model, our models can distinguish class I from class II MHCs and correlate their structural parameters with the TCR/pMHC’s potency to trigger T cell activation. The models are tested by mutagenesis using an MHC and a TCR mutated to alter conformation changes. The extensive comparisons between theory and experiment provide model validation and testable hypothesis regarding specific conformational changes that control bond profiles, thereby suggesting structural mechanisms for the inner workings of the TCR mechanosensing machinery and plausible explanations of why and how force may amplify TCR signaling and antigen discrimination.

its bound state during dissociation. Subtraction of this work from the interaction free energy tilts the energy landscape that governs the off-rate (Eq. 1, Fig. 1b). To calculate molecular stretch, we assume the TCR-pMHC complex to behave as a system of semi-rigid bodies of globular domains connected by semi-flexible polymers (Fig. 2b). As such, the total length change includes three components: 1) extension of individual globular domains, 2) various domain rotations about hinges, and 3) unfolding of secondary structures at specific regions.
For a globular domain without unfolding, the force-extension relationship is described by the three-dimensional freely-jointed chain model 1 : where i is end-to-end distance of the th-domain, B is the Boltzmann constant, T is absolute temperature, is force, and d ~ 100 pN is the elastic modulus of the folded globular domain 2 .
To calculate the work l ( ) , we project the above domain extensions onto the direction of force, which is taken as the z direction ( Fig. 2b) using two angles between the z axis and: 1) the normal direction of the bonding interface (θ) and 2) the line connecting the Cand N-termini of the MHC-I α3 domain excluding any unfolded residues (ϕ): Note that θ is a model parameter as it describes the tilting of the bonding interface as part of the force-induced conformational change, whereas ϕ is a model constant is measured from the crystal structure (Supplementary Table 1, see A.2 below).
We assume that partial unfolding in the molecular complex may occur at connecting regions of globular domains, in particular, the α1α2-α3 joint of the MHC-I and the Vα-Cα interdomain joint of the TCR (Fig. 2a). The former may be caused by dissociation of the noncovalent α1α2-β2m interdomain bond, which shifts the mechanical load originally borne by this bond to the α1α2-α3 hinge, resulting in its partial unfolding, as observed in SMD simulations 4 (Fig. 2a). Similarly, α1α2-β2m dissociation results in tilting of the bonding interface and load shifting from the Vβ-Cβ joint to the Vα-Cα joint, leading to partial unfolding of the latter joint (Fig. 2a). The unfolded polypeptides are flexible and can bear only tension but not moment, ensuring that their extension is along the direction of force, i.e., the z axis.
The force-extension relationship of the unfolded polypeptides can be described by an extensible worm-like chain (eWLC) model 5 : where p is the extension of the unfolded coil under force with the subscript indicating unstructured polypeptide, c = 0.36 nm and p = 0.39 nm are the average contour length and persistence length per unfolded amino acid, respectively 6-8 , p~ 50 µN is the elastic modulus of polypeptides 8 , j are polynomial coefficients for the improved approximation, and is the number of constituent amino acids in the unfolded polypeptide. In particular, we denote the respective numbers of amino acids in the unfolded MHC-I α1α2-α3 and TCR Vα-Cα joints to be p,MHC and p,TCR . Eq. 7 defines p as a function of f, which can be solved numerically to express in an explicit form: p c ⁄ = u,p ( ) = the extension per unit contour length for the polypeptide under force f.
Thus, the length of the TCR-pMHC complex at the transition state is ( Fig. 2b): Since we do not have prior knowledge about either number of unfolded amino acids, we will evaluate their sum * = p,TCR + p,MHC from curve-fitting of our model to the experimental data (see below). Since N ( ) is the length of the TCR-pMHC complex at the bound state ( Fig.   2b), we have ( ; 0) = N ( ). Finally, the integrand on the right-hand side of Eq. 3 can be written as

A.2. Simplifying assumptions and reducing parameters
The purpose of this subsection is to provide the details omitted in the main text of how  Table 1). This leaves only two structural parameters ( α3 , ) in our model to be evaluated by curve-fitting to data. Noting that the structure-based force function ( ; l ) scales with the characteristic extension change per unit change of molecular length such that l ( ) = ∫ z ( ) 0 , the dissociation rate of TCR-pMHC-Ⅰ bond can be written as: where z ( ) is given by Eq. 9.
To further reduce model parameters, we make additional assumptions as discussed in the main text. It is well-known that the fractions of free-energy change in biological interactions in liquid, such as unfolding and refolding of proteins, unbinding and rebinding of receptorligand bonds, and unzipping and rezipping of RNA or DNA, are small because of their limited dynamic transition time 9-12 . Such a rate limit results from the nature of biological interactions, e.g., polar/non-polar interactions, hydrophobic interactions, and charged interactions, which typically yield finite range of transition kinetics. This enables to roughly estimate the freeenergy barrier as ∆ 0 * ~l n( w / 0 ) where w~1 0 6 s -1 known as the prefactor 9-12 .

A.3. Defining the dissociation coordinate
The purpose of this subsection is to derive an operational way to determine the dissociation coordinate variable l used in the main text. We note that the total number of unfolded amino acids * is zero at the bound state before unfolding occurs, increases monotonically during progressive unfolding along the dissociation path, and reaches maximum at the dissociation point. Because * is not known a priori, it must be treated as a fitting parameter similar to α3 , , and 0 * . Since l is the contour length change along the dissociation path, we wish that l approaches its upper bound 0 * and depends on force implicitly through the model parameters * , α3 , and , as given below: where ∆( α3 , ) is the difference of the contour lengths except for the partially unfolded regions. Thus, Eq. 11 provides a constraint for l instead of introducing another model parameter. α3 and are determined for each * during the model fitting that searches for parameters to enable l → 0 * . Small errors, which may occur for the various contributions to ∆ (red lines in Fig. 2a labeled as force transmission lines), can be identified using a pair of ( α3 , ) values and the crystal structure for each complex. Specifically, average differences 〈∆〉 are -0.8 nm for the strong catch bond (SC) group ( α3 > 3 nm, > 30°), -0.4 nm for the weak catch bond (WC) group (1 nm < α3 < 3 nm, 8°< < 30°), and 0 nm for slip-only bond (SO) group ( α3 < 1 nm, < 8°), respectively (see Fig. 4e-i and associated text for the definitions of SC, WC, and SO groups). It has been well established that the contour length of a single amino acid is ~0.4 ± 0.02 nm/a.a 7,13,14 , implying that the model has a resolution of 2 amino acids. We further note that, even without conformational change, it is possible for the slip-only bond group to have 3 unfolded amino acids due to the limited resolution of the model.
Finally, the best-fit model parameters can be determined by finding the subset of best-fit parameters among possible * values that match the contour length change l to the freeenergy well width at zero-force 0 * , i.e., finding * such that l ( * , α3 , ) → 0 * (see Supplementary Table 2). Under this condition, ∫ z ( ) 0 = l ( ) → 0 * ( ) and Eq. 10 becomes identical to Eq. 2.

A.4. Model applications, curve-fitting strategies, and biological relevance
The purpose of this subsection is to outline the procedures of applying our model to The free-energy landscape can be constructed by substituting the best-fit model parameters into the following equations 14 : Thus, by using model parameters, the dissociation state coordinates relative to the bound state coordinates in the free-energy vs dissociation coordinate space can be defined as functions of force. Note that this force-induced change of the barrier height should be under the condition of small perturbation such that |− 0 * ( )| < ∆ 0 * . Since our fitting results show that the average free-energy barrier height at zero force is ~12 kBT ( 〈∆ 0 * 〉 ), the force range corresponding to a change of the barrier height of <10 kBT is reasonable for each dataset, i.e., force range corresponding to energy barrier heights in the range of 2 kBT < ∆ 0 * < 22 kBT.
Firstly, to fit the model-predicted reciprocal off-rate (Eq. 10) to the experimental bond lifetime vs force data, all four parameters were changed simultaneously to search for the minimum of the chi squared error. Second, every fitting curve was then determined by the parameter set that has the closest l value to 0 * (described as Eq. 11 in section A.3 and Supplementary Table 2). The fitting uncertainty of the best-fit parameters (or their ranges) were calculated by the differences between the previously obtained parameter set and parameter set from fitting to Mean ± SEM of bond lifetimes with the lowest RSS and the lowest Chi-square.
Of the 55 datasets analyzed, only one has 4 data points and this dataset shows slip bond; as such, is governed by two fitting parameters because the other two parameters are nearly zero.
All other datasets have 6-10 data points; therefore, over-fitting is not a problem.
To elucidate the biological relevance of the model parameters , * , Cαβ , 0 * , and ∆ 0 * , we examine their changes with varying bond lifetime vs force data obtained from different TCR and pMHC interactions that induce a wide range of biological responses. Finding correspondence between a group of model parameters individually and/or collectively with the biological response would be considered as support for the biological relevance of the model, because such correspondence suggests that the model can discriminate different TCR-pMHC interactions that result in differential T cell functions.

A.5. Class I model constraints
The purpose of this subsection is to check whether the model parameters obtained from data fitting are consistent with the constraints to which the model is subjected. Our experiments applied tensile force through the two ends of the TCR-pMHC complex, such that the force direction would always align with the line connecting the C-termini of the respective TCR and MHC molecules during dissociation, giving rise to the so-called pulling constraint. To formulate this pulling constraint in our model, we note that the pulling line is maintained so that the coordinate perpendicular to the force direction is invariant. As depicted in Supplementary Fig. 4a, several angle and length variables can be related using model parameters and structural constants by: where is the angle in an isosceles triangle constructed by rotating the α3 domain. By assuming that the α3 domain would be aligned with force, we estimate that the angle in the last Under this small angle assumption, Eq. 13a can be approximated by: Upon inversing Eq. 13b, we found that the tilting angle is a function of the end-to-end distance of the α3 domain, i.e., = ( α3 ). Setting = 25°, which seems reasonable as it approximates the maximum value of , the structural parameters obtained by fitting are scattered in-between two black curves on the α3 − plane marked as pulling constraint:  Fig. 4b). Thus, p,TCR can be simply calculated as p,TCR = * − p,MHC . Upon combining all information, the tilting angle ( TCR ) of variable domains of TCR can be described by simple trigonometrical function: where TCR is width of two interdomain hinges of the TCR α-and β-subunits measured from the crystal structures ( Supplementary Fig. 4a, tilting constraint). In this work we use TCR = 3.7 ± 0.3 nm as a representative width due to structure-to-structure variations. Thus, by comparing the tilting angle of the bonding interface (model parameter) to the tilting angle of the TCR TCR (derived from another model parameter p,TCR and structural constants TCR and c ), we can check the validity of tilting constraint using linear regression in the vs TCR plot ( Supplementary Fig. 4d).

B.1. Development, validation, and characterization
The purpose of this subsection is to present details of the development, validation, and characterization of the TCR-pMHC-II catch bond model omitted in the main text for simplicity, in a similar fashion as the TCR-pMHC-I catch bond model described in Section A.
The two models share exact the same framework but have different detailed form of the characteristic extension change z ( ) (Eq. 9). Comparing to the TCR-pMHC-I complex, the TCR-pMHC-II complex has different docking domains and pulling geometries (one vs two transmembrane domains on both TCR and MHC). For this reason, we assume that the forceinduced bonding interface tilting angle ( ) would be much smaller in the TCR-pMHC-Ⅱ than TCR-pMHC-I complex. The extension at the bound state can be defined as the end-to-end distance between both end-points identified by crystal structures (E8: 2IAM, 2IAN and 2B4 (as the substitution of 3.L2): 6BGA, 3QIB): where N is set to be 12.3 nm based on the crystal structures and linker ≈ 9.4 nm represents the linker (e.g., a leucine zipper) engineered at the C-termini of soluble pMHC-II constructs to stabilize both the MHC α-and β-subunits, which is often used in experiments for measuring where ∆( B.I , ) < 0.2 nm because < 10°. The best-fit model parameters can be determined by finding a subset of best-fit parameters among possible p,MHC values that match the contour length-change along the dissociation coordinate ( l ) to the width of the free-energy well at zero force ( 0 * ), i.e., finding p,MHC such that l � p,MHC � → 0 * (see Supplementary Table 5).
To characterize the class II model, we Fitting of class II model to experimental data and examination of the biological relevance of the best-fit parameters were done the same way as the class I model.

B.2. Class II model constraints
The purpose of this subsection is to describe the constraints of the class II model, which share similar ideas to those of the class I model (e.g., pulling and tilting constraints) but differ in their specific expressions. To formulate the pulling constraint, we again used the fact that the pulling force direction must aligns with the line connecting to the C-termini of the TCR and pMHC molecules so that the coordinate perpendicular to the force direction is invariant.
Using model parameters and structural constants, this pulling constraint can be written as: which can be solved for p,TCR explicitly: which can be approximately calculated using contour lengths of length components at forcefree state.
By assuming small angle perturbation, which seems reasonable, we further reduce the Thus, all terms including can be re-expressed by using Eq. 21.
Additionally, the tilting constraint can be expressed as follows: Thus, by using Eq. 19 and 21, only 4 fitting parameters, two structural parameters ( B.I , ) and two biophysical parameters ( 0 , 0 * ), were used to fit the class II model to data.

C. A general biophysical limit of model parameters
The purpose of this section is to describe a general biophysical limit that constrains the fitted model parameters, which is used in Fig. 6 to accept the correct model application to data of matched MHC class and reject incorrect model application to data of mismatched MHC class. The idea is that, even if the model is capable of fitting experimental data and the parameters are self-consistent with one another within the model, their values should be within known limits. A prototypical example of such a biophysical limit involves the molecular extension per unfolded amino acids. It follows from Eqs. 11 and 17 that the average molecular extension at zero force over all data (〈 0 〉) should be a linear function of * such that 〈 0 〉 = * + . The y-axis intercept can be determined from the slip bond data because slip bonds are not expected to have unfolded amino acids (i.e., * = 0) but still have a nonzero extension.
The slope is constrained by the fact that the contour length of a single amino acid has a small range (~0.4 ± 0.02 nm/a.a.) 7,13,14 ; thus, the average molecular extension (〈 0 〉) per unfolded amino acid should be bounded by: Imposing this range limit of 〈 0 〉 /a.a. would enable us to rule out inappropriate application of the model even if such application could achieve reasonable level of goodnessof-fit. Indeed, all results were below the biophysical limitation.
Conversely, a nearly zero estimate of would indicate that the model is inappropriate for catch-slip bond data because, for the model to fit catch-slip bonds, it requires * > 0 (see 2c-g and 3e, 3 rd row). A parameter estimation of ≈ 0 indicates the lack of dependence of the model behavior on * , which abolishes the model's ability to predict TCR signaling and antigen discrimination, making the model irrelevant to biology.
The sturdier class II than class I pMHC structure also precludes large rotation during conformational change at transition state, leaving only unfolding and stretching along the force.
Thus, the average molecular extension per amino acid should be close to 0.4 nm/a.a. On the other hand, ≈ 0 is expected from fitting the slip-only data because such data correspond to the * = 0 case, which makes it difficult to robustly estimate the correct value. Thus, the average molecular extension per amino acid should be well-correlated with each other (i.e., high level of goodness-of-fit as measured by 2 ) and in the range between 0 to 0.4 nm/a.a. In  Purified recombinant TCRs coated on beads were used in BFP measurement. * OT1 or 2C TCR expressed on CD4 + CD8 + thymocytes were used. † OT1, P14 or 2C TCR expressed on naïve T cells from transgenic mice were used. In other cases, 2C or 1G4 TCR expressed on 58α -1 β -1 hybridomas and J76 Jurkat cells were used. ¶ Soluble mouse N15 TCRαβ was used to measure bond lifetime in optical tweezers 19 . # Mouse TRBV TCRs (B13.C1 and B17.C1 with canonical docking orientation and B17.R1 and B17.R2 with reverse docking orientation) expressed on 5KC hybridomas interacted with NP366 bound to the D227K MT of H-2D b to prevent CD8 binding 18 ⁑ I1C mutation in β2m domain, G120C and C121S mutations in α domain (H2-K b ) were introduced. ** The recombinant TCR coated to the BFP bead was tested against pMHC expressed on a THP-1 cell. N.A. indicates not applicable.
Note that the best-fit parameters were extracted from equation (4) (two-pathway model) in the main text. All errors of fitting parameters are from the standard errors of bond lifetimes after optimizing the best-fit parameters.